#include "sphere.hpp"
namespace chen {

bool Sphere::interact(const Ray & _ray, SurfInteract * si) const {
    Vector3 oerr, derr;
    Ray ray = (*(this->w2o))(_ray, &oerr, &derr);
    //|o+td|^2 = r^2
    EFloat eox(ray.o.x, oerr.x), eoy(ray.o.y, oerr.y), eoz(ray.o.z, oerr.z), er(r);
    EFloat edx(ray.d.x, derr.x), edy(ray.d.y, derr.y), edz(ray.d.z, derr.z);
    EFloat t1, t2;
    bool sol = quadratic(
        edx*edx+edy*edy+edz*edz,
        (edx*eox+edy*eoy+edz*eoz)*(EFloat)2,
        (eox*eox+eoy*eoy)+(eoz*eoz-er*er),
        &t1, &t2
    );
    if(!sol) return false;
    if(t2.lower() <= 0 || t1.upper() > ray.t) return false;
    if(t1.lower() <= 0) {
        if(t2.upper() > ray.t) return false;
        t1 = t2;
    }
    _ray.t = t1.value;
    Point3 phit = ray(t1.value).resized(r);
    si->p = phit;
    Point3 pnorm = phit.normalized();
    si->ng = (Normal)pnorm;
    si->wh = ray.d.normalized();
    si->flipped = ray.o.lengthSquared() < r*r;
    Float l = std::sqrt(pnorm.x*pnorm.x+pnorm.y*pnorm.y);
    if(l == 0) {
        si->ngy = Normal(0, 1, 0);
    } else {
        Float frac = pnorm.z/l;
        si->ngy = Normal(-pnorm.x*frac, -pnorm.y*frac, l);
    }
    // Calculate differentials and surface local attributes
    Float cu = clamp(0, float_one_minus, (atan2(pnorm.x, pnorm.y)+(Float)(M_PI))/(Float)(2*M_PI));
    Float cv = clamp(0, float_one_minus, (asin(si->ng.z)+(Float)M_PI_2)/(Float)(M_PI));
    si->uv = {cu, cv};
    si->dpdu.z = 0;
    si->dpdv.z = (Float)M_PI*r*l;//(Float)M_PI*r*sin((Float)M_PI*v);
    // TODO correct abs value / possibly incorrect values
    Float dt0 = (Float)M_PI*r*pnorm.z;//(Float)(2*M_PI)*r*cos((Float)M_PI*v);
    Float dt1 = (Float)M_PI*r*l;//(Float)M_PI*r*sin((Float)M_PI*v);
    Float dt2 = -pnorm.y;//cos((Float)(2*M_PI)*u);
    Float dt3 = -pnorm.x;//sin((Float)(2*M_PI)*u);
    // TODO calculate dnduv
    si->dndu = si->dndv = {(Float)0, (Float)0, (Float)0};
    si->dpdu.x = -dt0*dt2;
    si->dpdv.x = dt1*dt3;
    si->dpdu.y = dt0*dt3;
    si->dpdv.y = dt1*dt2;
    // Here, be very cautious
    // If we use (const Shape*) to transform the this pointer
    // the virtual table offset will be lost, causing an unwanted
    // segmentation error in further code.
    // Is it a bug of gnu g++? 
    si->hit = (const Object*)this;
    si->p_err = abs((Vector3)phit)*gamma(5);
    *si = (this->w2o)->reversed(*si);
    return true;
}
bool Sphere::mask(const Ray & _ray) const {
    Vector3 oerr, derr;
    Ray ray = (*(this->w2o))(_ray, &oerr, &derr);
    //|o+td|^2 = r^2
    EFloat eox(ray.o.x, oerr.x), eoy(ray.o.y, oerr.y), eoz(ray.o.z, oerr.z), er(r);
    EFloat edx(ray.d.x, derr.x), edy(ray.d.y, derr.y), edz(ray.d.z, derr.z);
    EFloat t1, t2;
    bool v = quadratic(
        edx*edx+edy*edy+edz*edz,
        (edx*eox+edy*eoy+edz*eoz)*(EFloat)2,
        eox*eox+eoy*eoy+eoz*eoz-er*er,
        &t1, &t2
    );
    if(!v) return false;
    if(t2.lower() <= 0 || t1.upper() > ray.t) return false;
    if(t1.lower() <= 0) {
        if(t2.upper() > ray.t) return false;
        t1 = t2;
    }
    _ray.t = t1.value;
    return true;
} 
void Sphere::uniformSample (const Point2 & u, Interact * in, Point2 * uv, Float * pdf) const {
    Float x = cos(u.x*(Float)(M_PI*2.0)), y = sin(u.x*(Float)(M_PI*2.0));
    Float z = u.y*2-1;
    Float xyfrac = std::sqrt(max((Float)0.0, (Float)1.0-z*z));
    Point3 pnorm(x*xyfrac, y*xyfrac, z);
    Point3 hp = pnorm.resized(r);
    Vector3 p_err = (Vector3)abs(hp)*gamma(5);
    if(uv != nullptr) {
        Float cu = clamp(0, float_one_minus, (atan2(pnorm.x, pnorm.y)+(Float)(M_PI))/(Float)(2*M_PI));
        Float cv = clamp(0, float_one_minus, (asin(pnorm.z)+(Float)M_PI_2)/(Float)(M_PI));
        *uv = {cu, cv};
    }
    *in = w2o->reversed(
        Interact(
            hp,
            p_err,
            Vector3(0),
            (Normal)pnorm,
            this,
            false
        )
    );
    *pdf = (Float)1.0/area();
}

void Sphere::uniformSampleFrom (const Point2 & u, const Interact & from, Interact * in, Point2 * uv, Float * pdf) const {
    Point3 pfrom = (*w2o)(from.p);
    Float lsqr = pfrom.lengthSquared();
    if(lsqr < r*r) {
        uniformSample(u, in, uv, pdf);
        in->flipped = true;
        return ;
    } else {
        Float l = std::sqrt(lsqr);
        Normal nr((pfrom).normalized());
        Vector3 ny = nr.getPerpendicularUnit();
        Vector3 nz = cross(ny, (Vector3)nr);
        Float t = r*r/l;
        Float cosxtheta = (u.x*(r-t)+t)/r;
        Float sinxtheta = std::sqrt(max((Float)0, (Float)1-cosxtheta*cosxtheta));
        ny.resize(sinxtheta*r);
        nz.resize(sinxtheta*r);
        Vector3 nx = nr.resized(r*cosxtheta);
        Float xphi = u.y*(Float)(M_PI*2);
        Point3 pnorm(nx+ny*cos(xphi)+nz*sin(xphi));
        Point3 hp = pnorm.resized(r);
        *in = w2o->reversed(
            Interact(
                hp,
                (Vector3)abs(hp)*gamma(5),
                Vector3(0),
                (Normal)hp.normalized(),
                this,
                false
            )
        );
        if(uv != nullptr) {
            Float cu = clamp(0, float_one_minus, (atan2(pnorm.x, pnorm.y)+(Float)(M_PI))/(Float)(2*M_PI));
            Float cv = clamp(0, float_one_minus, (asin(in->ng.z)+(Float)M_PI_2)/(Float)(M_PI));
            *uv = {cu, cv};
        }
        *pdf = (Float)1/(((Float)1.0-r/l)*r*r*(Float)(2*M_PI));
    }
}


}